Exact numerical solution for a time-dependent fiber-bundle model 

with continuous damage 



O 

o 

Oh- 



L. Moral 

Departamento de Matemdtica Aplicada, Universidad de Zaragoza, 
50009 Zaragoza, Spain. 

J.B. Gomez 

Departamento de Ciencias de la Tierra, Universidad de Zaragoza, 
50009 Zaragoza, Spain. 

Y. Moreno 

The Abdus Salam International Centre for Theoretical Physics, 
Condensed Matter Group, P.O. Box 586, Trieste, 1-34014, Italy. 



o 



A.F. Pacheco 

Departamento de Fisica Teorica, Universidad de Zaragoza, 
50009 Zaragoza, Spain. 
(February 1, 2008) 



CO 



I 

C 

O 

o 



> 
o 

o 



C3 



-a 
c 

o 
o 



X 
S3 



A time-dependent global fiber-bundle model of fracture 
with continuous damage was recently formulated in terms of 
an autonomous differential system and numerically solved by 
applying a discrete probabilistic method. In this paper we 
provide a method to obtain the exact numerical solution for 
this problem. It is based on the introduction of successive 
integrating parameters which permits a robust inversion of 
the numerical integrations appearing in the problem. 

PACS number(s): 46.50.+a, 62.20.Fe, 62.20.Mk. 



I. INTRODUCTION 

Fracture in disordered media has for many years at- 
tracted much scientific and industrial interest Q-Q]. An 
important class of models of material failure is the fiber- 
bundle models (FBM) which have been extensively stud- 
ied during the past decades |f7|-[l^| . These models consist 
of a set of parallel fibers having statistically distributed 
strengths. The sample is loaded parallel to the fiber di- 
rection, and a fiber fails if the load acting on it exceeds 
a strength threshold value. When a fiber fails, its load 
is transferred to other surviving fibers in the bundle, ac- 
cording to a specific transfer rule. Among the possible 
options of load transfer, one simplification that makes 
the problem analytically tractable is the assumption of 
equal load sharing (ELS), or global load transfer, which 
means that after each fiber breaks, its stress is equally 
distributed among the intact fibers. Until very recently, 
the failure rule applied in standard FBM was discontinu- 
ous and irreversible, i.e., when the local load exceeds the 
failure threshold of a fiber, the fiber is removed from the 
calculation and is never restored. Recently, a novel con- 
tinuous damage law was incorporated into these models 



[OiLij. Thus, when the strength threshold of a fiber is 
exceeded, it yields, and the elastic modulus of the fiber 
is reduced by a factor a (0 < a < 1). Multiple yields of 
a given fiber are allowed, up to a maximum of n yield- 
ing events per fiber, where n is a small integer number 
which can be different for each fiber. This generalization 
of the standard FBM is suitable to describe a variety of 
elasto-plastic constitutive behaviors [ftE! -[l7| . 

The standard FBMs simulate the failure of a system at 
the microscopic level. Each fiber breakage can be mapped 
onto a new microcrack (with a typical size of a few /im) , 
or onto the extension of a previous microcrack. On the 
other hand, the continuous damage FBMs simulate fail- 
ure at a mesoscopic level. Now, each fiber in the model 
can be viewed as a small volume of the material. The 
term "small" depends on the size of the heterogeneities, 
but can be of the order of one millimeter for rocks. In 
each of these representative elementary volumes (REVs) 
in which the total volume can be divided, there are many 
potential sites for crack nucleation and growth, and the 
addition of each new crack will change continuously the 
elastic properties of the REV until its final failure when 
the accumulated damage surpasses a threshold. This 
threshold is identified in our model with the parameter n. 
Another important parameter in the model, the stiffness 
reduction factor, a, controls the amount of weakening 
that each yielding event introduces in a REV. The value 
a = 1 means no weakening, so that the elastic modulus 
of the REV remains the same irrespective of the number 
of yicldings, a rather unphysical situation. At the other 
extreme, the value a = means complete weakening after 
the first yield event. Thus, < a < 1 is the physically 
meaningful range for the stiffness reduction factor. In 
all the results given in the following sections, we have 
assumed that the initial elastic module of all the REVs 
is unity and that n is the same for all the fibers. The 
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randomness is incorporated in the REVs' liftimes, not in 
the elastic moduli. 

FBMs come in two settings, static and time-dependent 
or dynamic, and both of them have been applied to the 
standard and continuous damage settings [T^Jl^ , [l8| , [T9t . 
The static version of FBM simulates the failure of ma- 
terials by quasistatic loading. Drawing an analogy with 
what is carried out in a deformation experiment in the 
laboratory, a static FBM simulates a uniaxial or triax- 
ial, compressive or tensile, deformation test where the 
duration of the test is measured in seconds or minutes. 
In these models, the stress on each fiber is the indepen- 
dent variable and the strength of each element is the 
distributed random variable. On the other hand, the 
dynamic FBM simulates failure by creep rupture, static 
fatigue, or delayed rupture, i.e., a (usually) constant load 
is imposed on the system and the elements break because 
of fatigue after a period of time. The time elapsed until 
the system collapses is the lifetime of the bundle. Time 
acts as an independent variable, and the initial lifetime 
of each element, for a prescribed initial stress, is the inde- 
pendent identically distributed random quantity. Again, 
we can draw a clear analogy with a particular type of 
deformation experiments in the laboratory, the so-called 
creep experiments, where a heterogeneous material (rock, 
concrete, composite, ceramic alloy) is subjected to a con- 
stant or cyclic load, breaking after a period of time. The 
duration of these tests depends on the load imposed on 
the material and, more exactly, on the load compared 
to the short-term strength of the material (i.e., the load 
that causes the "instantaneous" failure of the same ma- 
terial in a fast uniaxial experiment). This load is usually 
expressed as a percentage of the short-term strength and 
the duration of the experiments is critically dependent 
on it. For rock, say, a sample will fail by creep after a 
few hours when subjected to a load 80% of the short- 
term strength, after a few weeks for a load 70% of the 
short-term strength, and after a few months or even years 
for lower working loads. The mechanism behind creep 
failure is subcritical crack growth, i.e., the slow exten- 
sion of microcracks with lengths smaller than the crit- 
ical crack length for instantaneous failure. Subcritical 
crack growth is due to a variety of processes operating 
near crack tips, the most important of them being stress 
corrosion, a chemical interaction between the crack tip 
and the environmental species, notably water, filling the 
microcracks that provokes the hydrolytic weakening of 
the atomic bonds of the material in the crack tip, where 
stress concentrations are highest. The crack propagation 
velocity is extremely sensitive to the applied load, sug- 
gesting exponential or power-law velocity functions with 
large coefficients or exponents. 

Indeed, in the dynamic FBMs the most widely used 
breaking rate function is the power law fl(i|-|l2|| , in which 
elements break at a rate proportional to a power of their 
stress, a p , where the exponent p is an integer called the 
stress corrosion exponent, for obvious reasons. This type 
of breaking rate will be assumed here and is another pa- 



rameter of the model. 

Our generalization of the dynamic global FBM |H| was 
restricted to the global transfer modality, and there we 
assumed that the size of the bundle, N , was very large. 
This enabled us to formulate the evolution of the system 
in terms of continuous differential equations. This type 
of equation, similar to those appearing in radioactivity, 
was first used by Coleman ||, and later in pjj| . In 18 we 
supposed an ELS bundle formed by N fibers which breaks 
because of stress corrosion under the action of an external 
constant load F ~ N ■ ctq, with <7o = 1. The breaking 
rate of the fibers, T, is assumed to be of the power-law 
type, T = <j p , f denotes the strain of the bundle and Y = 
1 represents the initial stiffness of the individual fibers. 
The original dynamic FBM was generalized by allowing 
one fiber to fail more than once, and thus we define the 
integer n as the maximum value of failures allowed per 
fiber. Besides, as mentioned before, the parameter a (< 
1) represents the factor of reduction in the stiffness of the 
fibers when they fail. As up to n partial yielding events 
are permitted per fiber, at any one time the population 
of fibers will be sorted in n + 2 lists. Thus N ~ Nq + 
Ni + . . . + N n + N', where N (i = 0,...,n) denotes the 
number of elements that have failed i times. TV' denotes 
the number of elements that have failed n + 1 times and 
therefore are inactive (i.e., they no longer support any 
load anymore). At t = 0, the N elements of the bundle 
form the list 0, N = N, and at t = T, N' = N. The 
specification, at a given time t, of the value of TVj, for 
i = 0, 1, ... ,n, provides the state of the system. In our 
continuous formulation the Ni are real positive numbers 
lower than N. 

As the external load F — N is supported by the present 
active fibers, we have N = f ■ (N + aNi + a 2 N 2 + . . . + 
a n N n ), and hence 



/ = N/(N + aN x + a 2 N 2 + ... + a n N n ). 
The time evolution equations are |lS|| : 
dN 



(1.1) 



dt 
dNt 

~df 

dN 2 

dt 



dNn 

dt 



f p (-N ), 
f{N - KNi), 
f"K{N x - kN 2 ), 

f P K n ^ 1 (N n -i — K,N n ), 



(1.2) 



where the ubiquitous constant factor n represents n = a p . 
This is an autonomous differential system. Its solution 
must fulfill the initial condition 



N (t = 0) = N 

JV 3 -(* = 0) = 0, j^0. 



(1.3) 



An alternative way of introducing a time-dependent 
rheological response in FBM is that of Cruz-Hidalgo et. 
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al. These authors incorporate a viscoelastic con- 

stitutive behaviour in their model through the mapping 
of each fiber to a Kelvin- Voigt element. They express 
the time evolution of the strain in each fiber by way of 
a differential equation. In their model, fibers break ir- 
reversibly when they surpass a statistically distributed 
strain threshold, whereas in our model multiple failures 
(yields) of a fiber are allowed, the variable which is statis- 
tically distributed is the lifetime of the fibers, and there 
is no explicit threshold dynamics. This different formu- 
lation implies that we can formulate the evolution of the 
system in terms of coupled differential equations, while 
the authors in ref. |19| have necessarily to use Monte 
Carlo simulations due to the lack of a global differential 
equation for the system. 



with 



In reference |18|], Eqs. (1.2) were solved by applying a 
numerical probabilistic method. The purpose of this pa- 
per i s to present an exact numerical method tha t solves 
Eq. (|L2l) , fulfilling the initial conditions (|L3|) . This 
method is explained in Section |n| In Section [01] wc 
present a discussion of the method and of the results. 
The reader will find a longer discussion of the physical 
results in Ref. |18|. This paper concentrates on the solu- 
tion method. 



II. EXACT NUMERICAL METHOD 

To simplify the notation, we first normalize the vari- 
ables 



Ni 

Xi = —, 1 = 0,1, 



(2.1) 



In terms of the Xi , the differential system to be solved is 
Xo = - f p x 



X j — f^f%? i^Xj — 'Y KiXj'j 

x (p) = l, x j (0) = 0, j = 1,2, 



(2.2) 



A dot on a variable means derivation with respect to 
time, and / and k are the same objects as in Section || 



!// = £' 



(2.3) 



i=0 



The system ( |2.2j ) admits a reduction of degrees of free- 
dom by eliminating t from the last n equations and by 
integrating with respect to xq: 



b (0) - 1 



6 w = 



In consequence, 



3-1 
1=0 



K 



j = 1,2, •••,71. 



(2.6) 



E l= ELo a i X 



fo(xo) 



l=i 



(2.7) 



Then, the first equation in (2.4) provides the relation xq 
versus t 



t = 



dxo 



x [fo(xo)] P X 



=0 a i x 



-dx (2.8) 



which, in principle, solves the problem because it relates 
t to xq and hence to any other Xj. However, the integral 
fl2.8p is, in general, improper for xq — » because the inte- 
grand is 0(xq K ~ ) and therefore the numerical relation 
t vs. xq is problematic. Specifically: 

a) If pK n — 1 < this integral is proper, 

b) if pK n — 1 > the integral is improper. 

Due to t he fa ct that the convergence occurs iff pK n — 1 > 
— 1, Eq. d2.8|) is always convergent, because in our model 
of fracture pn n > 0. 

Let e £ (0, 1); due to the fact that xo decays from 1 to 
0, ther e exists a time value to > such that xo{to) = £■ 
If (2.8) is improper, we perform the following change of 
parameter: xq = yo — * y\, such that 



xo = -f p x 

ij = j^kP [Xj— i ~ 
with x (t ) = e, yi(t ) 



txj]; j = l,2,---,n, (2.9) 
1, and t > to- From here 



dxo 
dyi 

Hence 



KXo 

yi 



1/k 

xo = ciy^ , 



ci = (x (t )) = e. (2.10) 



xo = - fxo, 
dxj K,3~ 1 (nxj — Xj-i) 
dxo 



xo 



From Eq. (2A) we obtain 

i 

Xj = ^2 by'xfi , i = 0,1, 



(2.4) 



(2.5) 



f=-h(yi) 



1=0 



(2.11) 



1=0 



(2.12) 
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With P^' = &£V , n 

0,l,---,n). 

In these circumstances ( p.ll| ) and the equation 



(i) 



( 



t-to = 



dyi 



V1 (/i(yi)) yi 



2/1 



(2.13) 



describe t, xq, • • • ,x n in terms of the y\ parameter, for 

t>h. 

As the integrand of fl2.13|) is 0(y pK ~ 1 ), then 

a) If — 1 > (|2.13| ) is a proper integral, 

b) if pK n ^ 1 — 1 < ( |2.13 ) is an improper integral, but 



(2.13) is always convergent. 

Now, as y± decays to zero, there exists a time instant 
ti > t such that y\(ti) = e. And by considering the 
change of parameter yi — > 2/2 given by the conditions 



*0 = -f P XQ 

yi = -Kfyi 
2/2 = -k 2 / p 2/2 

xj = f p K? \Xj—\ ~ KX j] ! J = 1) 2, ■ • • , n, 
with 2/2(^2) = e, 2/2 (ii) = lj and t > t\, we have 



(2.14) 



*/2 



«2/i 

2/2 

and hence 



2/i = C22/2 



2 1A \ C2 = (!/!(*!))=£, (2.15) 



J=0 



(=0 



/ =: /2(2/2) 



(2.16) 



(2.17) 



with identical meaning as before for pf and a 
Then, ( |2.16| ) and 



(2) 



t-tl = 



dy2 



(2) K *-2\ p 
2/2 



2/2 



/ w (/2(2/2)) P 2/2 

(2.18) 

describe xq, • • ■ , x n m terms of 2/2 > for t >t\. Besides, 



as the integrand of (2.18) is 0(y% K ), then 



2.18| ) is a proper integral, 

is an improper integral, but 



a) if pn n - 2 - 1 > 

b) if pn n - 2 - 1< ( pJ8 
always convergent. 

The process followed so far is generalized in the way 
expressed in Table | where in the end 



/ =: fn{y n ) 



T n a {n) v K ' 



(2.19) 



and therefore 



t — in— 1 



c?2/« 



Vn {fn{y n )) p yn 



1 ( V!!' u n •//. 



I 071 



-dy n 
(2.20) 



whose integrand is 0(y? x ), that is, integral (2.20) is 
always proper. 



III. RESULTS AND CONCLUSIONS 

The simple formalism written in Section || can be 
expressed, for example, in a brief program of MATH- 
EMATICA and its results graphically appreciated. We 
omit here the program but it can be provided on request. 
By fixing the constants at the following values: n = 3, 
a = 0.6, p = 2, and e = 0.1, in Fig. 1 the value of the 
working parameters 2/i are represented vs. time. Note 
that their range of definition is from 1 to e, except for 2/3 
which ends at for t 3 = T, i.e., the actual lifetime of the 
bundle. 

In Fig. 2 we again show the evolution of the working 
parameters and also the evolution of the four lists Xi of 
elements in the problem. 

The strategy developed in Section [n| can be summa- 
rized in a few sentences. First, let us observe Fig. |^ to 
appreciate the time evolution of the different lists: while 
Xo monotonously declines from 1 at t = to at t = T, 
the lists xj, j — 1,2,3 start from at t = 0, rise to a 
maximum and then monotonously decline to at t = T 
(strictly speaking, all the lists vanish at the same time). 
The last list j = n is special in the sense that it is the 
only one that tends to with an infinite slope when t 
tends to T. 

The analytical resolution of Eq. ( |2.2|) is impossible be- 
cause of the nonlinearity introduced by the f p factors. 
This source of complexity is partly overcome after hav- 
ing recognized the par tial reduction of degrees of free- 
dom expressed in (2.4). This partial reduction leads to 



the re lati on between Xj, j — 1,2, •••,n and xq, hence 
from ( |2.8| ) one has solved in principle the time evolution 
of Xo, and of the rest of Xj. But, in (2.8) one also rec- 



ognizes that this integral is improper. This is the real 
problem we face for the numerical inversion t «-> xq in 
the region where Xn is very small. In intuitive terms, 
this shows in Fig. |2] because beyond a certain time, xo 
is no longer significant and its relation with t becomes 
"delicate". Therefore we have used xo — j/o a s a good 
integration parameter only up to t = to. Beyond this 
point we successively introduced other "artificial param- 
eters" 2/1, 2/2, • ■ • ! 2/n which in the corresponding time in- 
terval play the role performed by xq from to to- Using 
these parameters, we are able to robustly relate all the 
variables Xi to t in the whole interval from to T. 

At the end of the process, the last integral is always 
proper, which allows a robust numerical inversion in the 
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vicinity of t = T. Intuitively, this is clear in Fig. || where 
we appreciate the abrupt fall-off of X3. 

In th e com ments writte n in Section || after Eqs. ( |2.8| ), 
( 2.13 ), ( 2.18 ), and ( 2.20 ) regarding the nature of those 
integrands, we noted that in general they behave as 
0(yf K This implies that the condition 



pn 



tells us the value of i = i c , 



i r > n — 



lnp 
Un~d' 



(3.1) 



(3.2) 



such that, for i > i c , the respective integral is proper and 
there is no need to introduce more artificial integrating 
parameters. 

The reader should note that the e introduced in the 
method is not a limiting factor of precision, but merely 
sets the temporal ranges of the various integrating pa- 
rameters yi. In our procedure, the only source of inac- 
curacy is the precision of MATHEMATICA, used for the 
numerical inversion of the integrals. 

As a final conclusion we would say that the exact nu- 
merical method presented in this paper to solve this fiber- 
bundle problem does not predict any new qualitative re- 
sult with respect to what was obtained using the approx- 
imate method of Ref. |18|. Therefore, no new physical 
conclusions can be drawn from here. 

The use of this exact strategy in other scientific prob- 
lems that are cast as an autonomous differential system 
will be considered in the next future. In this respect, 
clear candidates are some ecological problems and mod- 
els of infection spreading fl2Q,M . 
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FIG. 1. Time evolution of the four integrating parameters 
J/o, y\ , J/2 and 1/3 for a system with n — 3, a = 0.6 and p — 2. 
Note that their range of definition is from 1 to e, except for 
1/3, which goes from 1 to 0. 
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FIG. 2. Time evolution for the four integrating parameters 
and the four variables xo, X\, %2 and x% for a system with the 
same parameters as for Fig. 1. 



TABLE I. General terms in the procedure. 



Time interval 



Condition 



Parameter 



[0,to] 
[toM 

[tiM 



3/o (to) = e 

yi(*i) = e 



2/o = x 

j/i such that yi = —nfyr, Vi(to) = 1 
j/2 such that t/2 = -K 2 f p y 2 ; yn(ti) = 1 



J/n-l(*n-l) = e 



j/„ such that y n = —K n f p y n ; y n (t n -i) = 1 



